#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <argtable2.h> //command line argument parcer
#include <Magick++.h> //imageMagick for c++

using namespace std;
using namespace Magick;

/* global arg_xxx structs */
//struct arg_lit *x,*y,*z;
struct arg_int *x,*y,*z,*argImageQuality;// *scal;
struct arg_int *argDimension;
struct arg_dbl *argResolution, *argOffset;
struct arg_file *o, *i;
struct arg_end *end;
struct arg_rem *R1;
struct arg_rem *argR0,*argR1,*argR2,*argR3;
struct arg_lit *HELP, *dataIsInt, *globalNorm, *argPNG, *argColor;

///MAKE PARALLEL WITH PTHREADS!!
///Slice Direction!!
///Transparency Value
///Set intensity/color scale

int main(int argc, char **argv)
{
/* the global arg_xxx structs are initialised within the argtable */
	void *argtable[]={
		HELP = arg_lit0(NULL, "help", "Show this help menu"),
		x = arg_int0("xX", NULL,"<nx>", "the x dimension of raw data"),
		y = arg_int0("yY", NULL,"<ny>", "the y dimension of raw data"),
		z = arg_int0("zZ", NULL,"<nz>", "the z dimension of raw data"),
		argImageQuality = arg_int0("qQ","ImageQuality","<n>","PNG: 100 for best compression."),
		R1 = arg_rem(NULL,"JPG: 100 for lossless. (Default: 75)"),
//		scal = arg_int0(NULL, "scalar","<n>", "foo value"),
		dataIsInt = arg_lit0("I", "int,Int", "Process integer data (default: float)"),
		globalNorm = arg_lit0("gG", "GlobalNorm", "Use constant intensity scaling (default: per-layer)"),
		argColor = arg_lit0("c","colorMap","output image as color map (jet)"),
		argPNG = arg_lit0("pP","PNG,png", "Lossless image compression. (Default: JPG)"),
		i = arg_file1("i","InFile","<file>", "binary array input file"),
		o = arg_file1("o", "OutFilePrefix","<file>", "output file(s) base name"),

		argR0 = arg_rem(NULL,"[Options below typically in .xargs files]"),
		argDimension = arg_intn(NULL, "Dimension","<x>",0,3,"Number of image data elements (voxels)."),
		argR1 = arg_rem(NULL,"* repeat argument for <y> and <z> in order."),
		argResolution = arg_dbln(NULL, "Resolution","<x>",0,3,"(ignored) Physical size of image voxels (cm)."),
		//argR2 = arg_rem(NULL,"* repeat argument for <y> and <z> in order."),
		argOffset = arg_dbln(NULL, "Offset","<x>",0,3,"(ignored) Image offset from isocenter (cm)."),
		//argR3 = arg_rem(NULL,"* repeat argument for <y> and <z> in order."),

		end = arg_end(20)
	};

	if (arg_nullcheck(argtable) != 0)
		printf("error: insufficient memory\n");

	int nerrors = arg_parse(argc,argv,argtable);

	//help message
	if( HELP->count ){
		printf("\nUsage: bin2jpg");
		arg_print_syntax(stdout,argtable,"\n\n");
		arg_print_glossary_gnu(stdout,argtable);

		arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
		return 0;
	}

	if (nerrors==0)
	{
		//print parced inputs
		//printf("-x = %d\n", x->ival[0]);
		//printf("-y = %d\n", y->ival[0]);
		//printf("-z = %d\n", z->ival[0]);
		//printf("-i = %s\n", i->filename[0]);
		//printf("-o = %s\n", o->basename[0]);

		int rx;
		int ry;
		int rz;
		//assign resolution variables
		if(x->count && y->count && z->count){
			rx = x->ival[0];
			ry = y->ival[0];
			rz = z->ival[0];
		}else if( argDimension->count == 3){
			rx = argDimension->ival[0];
			ry = argDimension->ival[1];
			rz = argDimension->ival[2];
		}else{
			printf("ERROR: Please specify geometry, see --help.\n");
			arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
			return 0;
		}
		int numberOfElements = rz*rx*ry;

		//open input file else error
		FILE *inputFile;
		inputFile = fopen(i->filename[0],"rb");
		if( inputFile == NULL )
		{
			printf("Can not open %s for input.\n",i->filename[0]);
			return 1;
		}

		//obtain file size
		fseek (inputFile , 0 , SEEK_END);
		long lSize = ftell(inputFile);
		rewind (inputFile);

		//check file size
		printf("%s is %li bytes\n",i->filename[0],lSize);
		if( lSize != (long)sizeof(float)*numberOfElements )
		{
			printf("File size is not consistant with input parameters.\n");
			return 1;
		}

		InitializeMagick(NULL);

		//allocate array in memory
		float* buffer = (float*) malloc(lSize);
		if (buffer == NULL) {fputs("Memory error",stderr); exit(2);}

		//handle integer input (simple)
		if( dataIsInt->count ){
			int* tempBuff = (int*) malloc(lSize);

			size_t result = fread (tempBuff,1,lSize,inputFile);
			if ((long)result != lSize) {fputs("Reading error",stderr); exit(3);}
			lSize/=sizeof(float); //now lSize is number of elements
			for(int v=0; v<lSize;v++){
				buffer[v] = float(tempBuff[v]);
			}
		}else{
			//read array
			size_t result = fread (buffer,1,lSize,inputFile);
			if ((long)result != lSize) {fputs("Reading error",stderr); exit(3);}
			fclose(inputFile);

			 if( globalNorm->count ){
				lSize/=sizeof(float);
			 }

		}

		// for better normalization :
		//find maximum value for each frame (z-layer)
		float* normFactor = (float*) malloc(sizeof(float)*rz);
		float minValue;// = (float*) malloc(sizeof(float)*rz);
		float* offsetValue = (float*) malloc(sizeof(float)*rz);

		//read in color map:
		//Image myMap;
		//myMap.read("jetHi.png");

		//for low-level color mapping
		Color myMap[10];
		float HueScale = 0.834f/5.0f;
		for(int i=5; i>=0; i--){
			myMap[i] = ColorHSL(0.166+HueScale*float(i),1.0,0.5);
		}


		if( dataIsInt->count || globalNorm->count ){ //using global normalization
			minValue = *min_element(buffer,buffer+lSize);
			float maxValue = *max_element(buffer,buffer+lSize);
			normFactor[0] = 1.0f / (maxValue - minValue);
			//avoid -0.0
			if(minValue == 0.0f){
				offsetValue[0] = 0.0f;
			}else{				
				offsetValue[0] = -minValue;//*normFactor[i];
			}
			printf("Scaling all images from 0 to %f\n", maxValue);
		}else{ //using per frame normalization
			for(int i=0;i<rz;i++)
			{
				minValue = *min_element(buffer+i*rx*ry,buffer+(i+1)*rx*ry);
				//maxValue = *max_element(buffer+i*rx*ry,buffer+(i+1)*rx*ry);
				normFactor[i] = 1.0f / (*max_element(buffer + i * rx * ry, buffer + (i + 1) * rx * ry) - minValue);
				//avoid -0.0
				if(minValue == 0.0f){
					offsetValue[i] = 0.0f;
				}else{				
					offsetValue[i] = -minValue;//*normFactor[i];
				}
				//printf("max found: %f\n", maxValues[i]);
			}
		}

		//save a jpeg file for each layer
		float* tmpPtr;
		char tempFileName[512];
			//(char*) malloc(sizeof(char)*strlen(o->basename[0])+(int)ceil(log10((float)rz))+1);
			//make file name array long enough to handle the added numbers
		//Image slice_image( Geometry( (unsigned int)rx, (unsigned int)ry ), Color(MaxRGB, MaxRGB, MaxRGB, 0));
		for(int frame=0; frame<rz; frame++)
		{
			tmpPtr = buffer+frame*rx*ry;

			if( dataIsInt->count || globalNorm->count ){ //for global normalization
				for(int i=0;i<rx*ry;i++)
				{
					tmpPtr[i] = (offsetValue[0] + tmpPtr[i])*normFactor[0];//try this on GPU? (FMA)
				}
			}else{ //using per frame normalization
				for(int i=0;i<rx*ry;i++)
				{
					tmpPtr[i] = (offsetValue[frame] + tmpPtr[i])*normFactor[frame];//try this on GPU? (FMA)
				}
			}

			//Blob blob(tmpPtr, rx*ry );
			Image slice_image(rx, ry, "I", FloatPixel, tmpPtr); //resolution, I for grey data, data type, array!
			//slice_image.size("512x512");
			//slice_image.magick( "I" );
			//slice_image.read(blob);
			
			
			//unsigned int constant: MaxRGB
			//apply color map
			if( argColor->count){
			
				//turn on alpha channel:
				slice_image.type(TrueColorMatteType);
				//slice_image.map(myMap,true);
				//Quantum R=MaxRGB,G=0,B=0,A=MaxRGB;
				int idx;
				//float maxFloat = float(MaxRGB);
				
				//try steps use HSV and hue starting at 85 to 256
				for(int i = 0; i<rx; i++){
					for(int j=0; j<ry; j++){
						idx = rx*j+i;
						//transparency threshold
						if(tmpPtr[idx] <= 0.2f){ ///TODO: make this an argument
							slice_image.pixelColor(i,j,Color(0,0,0,MaxRGB));
						}else if(tmpPtr[idx]<= 0.4f){
						 	slice_image.pixelColor(i,j,myMap[0]);
						}else if(tmpPtr[idx]<= 0.6f){
						 	slice_image.pixelColor(i,j,myMap[1]);
						}else if(tmpPtr[idx]<= 0.8f){
						 	slice_image.pixelColor(i,j,myMap[2]);
						}else if(tmpPtr[idx]<= 0.9f){
						 	slice_image.pixelColor(i,j,myMap[3]);
						}else if(tmpPtr[idx]<= 0.95f){
						 	slice_image.pixelColor(i,j,myMap[4]);
						}else{
						 	slice_image.pixelColor(i,j,myMap[5]);
						}
					}
				}
			}
			
			
			if( argImageQuality->count){ //set image quality
				slice_image.quality(argImageQuality->ival[0]); //100 is best quality, 75 is default
			}
			if( argPNG->count ){
				slice_image.magick("PNG");
				sprintf(tempFileName,"%s%i.png",o->filename[0],frame);//o->basename[0],
				slice_image.write(tempFileName);
				printf("\rSaved: %s", tempFileName);
				//printf("\rProgress: (%.0f%%)", float(frame+1)/float(rz)*100.0f);
				fflush(stdout);
			}else{
				slice_image.magick("JPEG");
				sprintf(tempFileName,"%s%i.jpg",o->filename[0],frame);//o->basename[0],
				slice_image.write(tempFileName);
				printf("\rSaved: %s", tempFileName);
				//printf("\rProgress: (%.0f%%)", float(frame+1)/float(rz)*100.0f);
				fflush(stdout);
			}
		}
		printf("\nDone!\n");
	}
	else
		arg_print_errors(stdout,end,"gDPM");


	arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));

	return 0;
};
